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Background 

Some  refinements  were  described  in  Ref.  1  for  the  ship  tracking  algorithm 
developed  in  Ref.  2.  The  availability  of  these  two  references  is  assumed  in 
this  report.  One  of  the  refinements  of  Ref.  1  is  extended  here  to  provide 
better  responsiveness  to  strong  target  maneuvers.  Also,  another  modification 
is  introduced  to  allow  the  tracking  algorithm  to  operate  on  bearing-only 
reports. 

Upon  the  receipt  of  each  new  position  report,  the  tracking  algorithm 
basically  operates  in  planar  (x,y)  coordinates  by  first  updating  an  estimate 
of  the  ship's  motion-state  vector  (x.y.x.yj  with  a  Kalman  filter,  then 
updating  an  estimate  of  a  2x2  "driving  noise"  covariance  matrix  Q  with  the 
"innovations"  from  this  filter.  This  Q  matrix  estimate  is  then  modified  to 
make  it  diagonal  with  respect  to  the  current  estimate  of  the  velocity  vector. 

The  main  refinement  of  Ref.  1  was  to  adjust  the  velocity  components  of  the 
"state  covariance  matrix"  at  this  point  whenever  the  estimate  of  the  driving 
noise  had  increased.  In  the  currently  estimated  "in-track"  and  "cro3s-track" 
coordinates,  this  adjustment  consisted  of  an  increase  in  the  (2x2)  velocity 

covariance  submatrix  by  (t-t  )_1Q,  where  each  component  of  Q  is  the  greater  of 

o 

zero  and  the  change  in  the  corresponding  component  of  Q  from  its  preceding 

value,  and  where  t  is  the  current  time  and  t  is  the  time  that  tracking 

o 

started.  This  refinement  produced  dramatic  performance  improvements  in  some 
cases  by  correcting  the  original  algorithm's  tendency  to  adhere  too  rigidly  to 
the  "average  velocity"  estimated  from  early  data.  The  rationale  for  this 
particular  adjustment  of  the  velocity  covariance  submatrix  was  that  it  is 
exactly  the  correction  that  would  be  needed  in  the  state  covariance  matrix 
after  two  position  reports  if  this  matrix  were  generated  by  the  Kalman  filter 
using  too  low  a  value  of  the  driving  noise. 

For  the  third  and  following  position  reports,  this  same  rationale  is 
extended  here  to  derive  corresponding  adjustments  of  the  (2x2)  position 
covariance  submatrix  and  the  expected  values  (Kalman  filter  estimates)  of  the 
position  and  velocity.  These  further  adjustments  do  not  make  the  kind  of 
dramatic  improvement  produced  by  adjusting  the  velocity  covariance  submatrix, 
but  they  do  cause  the  algorithm  to  respond  more  rapidly  and  accurately  to 
strong  target  maneuvers.  This  reduces  the  need  for  adopting  ad  hoc  maneuver 
detection  schemes,  which  often  have  unclear  effects  on  the  "error  ellipses" 
and  location  likelihoods  that  can  be  generated  by  this  algorithm. 
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Also  developed  here  is  an  alternate  way  to  updating  the  estimate  of  the 
driving  noise  matrix  with  the  Kalman  filter  residuals.  This  alternative  is 
based  on  the  approximation  that  the  in-track  and  cross-track  components  of  the 
driving  noise  are  equal  in  intensity  and  independent.  This  simplification  is 
essential  in  using  bearing-only  reports  (or  position  reports  with  long,  narrow 
"containment  ellipses")  to  estimate  the  driving  noise  because  the  procedure 
used  in  Ref.  1  depends  heavily  on  the  individual  reports  localizing  the  ship 
in  two  dimensions.  It  could  probably  also  be  used  with  well-localized  reports 
without  sacrificing  very  much  precision,  since  the  estimated  driving  noise  in 
such  cases  is  usually  fairly  isotropic  anyway  when  the  procedure  of  Ref.  1  is 
used. 

Only  planar  tracking  in  x,y  coordinates  is  discussed  here.  Extensions  to 
tracking  on  a  sphere  can  be  made  along  the  lines  described  in  Ref.  1  and  2. 

New  Method  for  Altering  Posterior  State  Density  Parameters  in  Response  to  an 


Increase  A  in  Estimated  Driving  Noise  Intensity  Matrix 


In  the  following,  let  M  denote  the  posterior  state  covariance  matrix  just 
before  the  report  at  time  t^  and  let  P  denote  this  matrix  just  after  the  pre¬ 
ceding  report.  Also  let  M  and  P  be  partitioned  into  2x2  partitions 
corresponding  to  position  and  velocity  variables  as 


M1  i  m2 

—f~r - 

w  ^  )  w 


and  P  s 


P1  i  p2 


2  !  ‘3 


where  M^  is  a  position  covariance  matrix,  etc.  If  (and  only  if)  a  po3itive- 

semidefinite  increase  A  occurs  in  the  processing  of  the  report  at  time  t^ 

(i.e.,  after  the  Kalman  filter  updating  of  the  posterior  state  density 
parameters) ,  these  parameters  are  then  further  altered  in  this  new  method  as 
follows: 

1.  Velocity  Covariance  Submatrix  K: 

For  the  second  and  following  reports,  K  is  replaced  by  K  +  5K,  where 

5fC  '  t  “  '  A’ 
i  o 


with  a  certain  positive  definite  upper  limit  on  sK.  t  is  the  time  of  the 

initial  report.  This  adjustment  is  the  refinement  of  Ref.  1  mentioned 
earlier,  and  is  described  in  detail  there. 


2.  Velocity  Mean  Q  (a  2-vector) : 
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For  the  third  and  following  reports,  it  is  assumed  that  also  should 
have  been  increased  by  SK  and  that  the  value  is  negligible.  Hence,  for  the 
third  and  following  observations  only,  v  is  replaced  by  v  +  5$,  where 

T  1 
0  +  5$  =  v  +  (M1  +  A  f-~t  )  <Mi  +  At  +  R)  (z  -  x), 

i”  o 

according  to  the  Kalman  filter  equations,  and  where  t  =  t  -t  and  the  rest 

l  i-1 

of  the  notation  is  as  in  Refs.  1  and  2.  For  simplicity  and  from  good  numeri¬ 
cal  experience,  the  approximation 

v'sv+M^  (M^  +  R  +  At)(z-x), 

is  adopted.  Hence,  the  correction  for  v  becomes 


3.  Position  Covariance  Submatrix  B: 


No  good  theory  i3  available  for  adjusting  the  entire  (4x4)  state  covar¬ 
iance  matrix  for  large  A,  which  is  the  case  of  major  interest,  because  of 
incomplete  state  observability.  Hence,  the  position-velocity  submatrix  is  not 
adjusted,  and  a  positive-semidefinite  (2x2)  increment  5B  is  sought  for  B,  so 
that  the  resulting  (4x4)  state  covariance  matrix  is  still  positive-definite. 
Now 

B  *  M  -  N  (M  ♦  R) -1M 1 . 

After  an  estimated  A  >  0,  we  only  consider  one  previous  report  because  posi¬ 
tion  is  directly  observable.  Hence,  from  the  Kalman  filter  propagation  equa¬ 
tions,  the  difference  caused  by  Q  being  augmented  by  A  is 

SM  ^  z  At  , 

SO 

B+$B  =  M^  +  At  -(M1  +  At)  (M  ^  ♦  R  +  At)"1(M1  +  At). 

If  terms  of  second  order  in  At  are  neglected  here, 
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5B  =  At  +  M1(M1+R)~1M1-(H1+At)[(M14.R)"1-(M14.R)"1At(M1-*.R)"1](M1+At) 
=  At  -  At (M  R)_1M  j-M  1  (M  ^R)'1  At  +  M  ^ (M  1+R) _1  At (M  1+R)-1M  1 
=  tI-M1(M1+R)-1]AT[I-(M1+R)'1M1] 

=  R(M1+R)'1At(M1+R)"1R 


Now,  for  one-dimensional  motion, 


CM  ”  +  R)(M  +  R  +  At) 


exactly  for  any  At-  Hence,  to  be  conservative  and  to  compensate  for  the  fact 
that  the  position-velocity  submatrix  is  not  being  altered  here,  we  use 


6B  =  R(M1+R+At)-1At(M1+R+At)-1R. 


This  adjustment  is  also  not  used  until  the  third  report  because  the  relatively 
large  initial  speed  variance  gives  good  results  without  it  at  the  second  up¬ 
dating  time. 


4.  Position  Mean  8:  (2-vector) 

From  the  Kalman  filter  equations  for  updating  with  the  position  report  z 
at  time  t^ , 

8  =  x  +  M  ^M  1+R)”1(z-x) 

and 

z  -  l  =  tl-M  (M1+H)”1](z-x>. 

As  before, 

«M1  =  At 

is  taken  after  an  estimated  positive-definite  A.  Thus  8  is  replaced  by 
8  ♦  4$,  with 


4 


Now, 


2  +  «2  =  x  +  (M  +At)(M  +At+R)_1(z-x). 


z-(2+«2)  =  [I-(M1-fAr)(M1+AT+R)"1](z-x) 

=  R  (M  +At+R) ’(z-x) 

=  R(M  ^At+R)"1  (M  ^At+R-At)  (M  1+R)_1  (z-x) 

=  RtI-(M1+AT+R)"1Ar](M1+R)_1(z-x) 

=  CM1+R-M1-R(M1+At+R)“1At](M1+R)“1(z-x) 

=  [I-M1(M1+R)"1](z-x)-R(M1+R+At)At(M1+R)“1(z-x) 

=  (z-x)-R(M1+R+At)At(Mi+R)~1(z-x) . 

Thus,  the  adjustment 

<52  =  R(M1+At+R)At(M1+R)"1(z-x) 

is  used,  but  not  until  the  third  report  is  processed  for  the  same  reason  as  in 
the  case  of  the  position  covariance  submatrix. 


Summary 


Only  adjustment  1  is  made  after  the  second  position  report  (and,  of 
course,  only  if  a  positive-definite  A  is  estimated  then  by  the  procedure  des¬ 
cribed  in  Ref.  1);  then  all  of  the  adjustments  1-4  are  used  thereafter  when¬ 
ever  A  >  0  is  estimated.  Making  these  adjustments  only  when  A  >  0  is  compen¬ 
sated  by  the  fact  that  corrections  are  only  made  one  step  back  for  the 
position-component  parameters  B  and  2  of  the  posterior  probability  density 
function. 

An  example  of  the  numerical  significance  of  using  all  four  of  these  cor¬ 
rections,  instead  of  only  the  first  one  as  was  done  in  Ref.  1,  is  shown  in 
Figures  1  and  2.  The  input  report  data  is  identical  in  both  cases,  and  cor- 
responds  to  a  ship  making  a  135°  turn  and  a  speed  change  at  time  t^.  To  keep 
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the  figures  simpler,  the  two-sigma  (90%  containment)  ellipses  for  the  output 
location  estimates  are  displayed  only  for  times  near  this  maneuver,  which  is 
when  they  are  most  different  for  the  two  tracking  algorithms.  The  further 
corrections  instituted  here  do  improve  the  tracking  performance  following  this 
maneuver,  especially  in  the  updated  containment  ellipse  at  the  sixth 
observation  time  and  the  projected  containment  ellipse  at  the  seventh 
observation  time.  The  containment  ellipses  are  also  somewhat  tighter  (but 
still  covering  the  target)  near  the  end  of  this  track  when  the  new  corrections 
are  added,  but  these  are  not  shown  in  the  figures. 


Alternate  Method  of  Estimating  Driving  Noise  Intensity 


The  basic  tracking  algorithm  here  treats  ship  motion  as  a  steady  average- 
velocity  drift  superimposed  on  a  two-dimensional  Brownian  motion  due  to 
maneuvering.  Beth  the  average  velocity  (a  2-vector)  and  the  intensity  para¬ 
meter  of  the  Brownian  motion  (the  2x2  "driving  noise  intensity"  matrix)  are 
estimated  from  the  incoming  position-report  data.  The  procedure  advocated  in 
Ref.  1  estimates  this  intensity  as  if  the  in-track  and  cross-track  components 
of  the  Brownian  motion  maneuvers  were  statistically  independent.  This  simpli¬ 
fication  was  adopted  because  such  independence  would  normally  apply  to  ship 
motion,  and  because  the  report  data  would  not  usually  support  a  usefully  reli¬ 
able  estimate  of  anything  more  refined  than  this.  These  two  manuevering  com¬ 
ponents  were  not  presumed  to  be  of  equal  intensity,  however,  because  it  was 
thought  that  allowing  their  estimates  to  be  different  would  often  lead  to 
significantly  more  precise  estimates  of  ship  motion ;  for  example,  when  a  ship 
zigzags  about  an  average  course,  or  when  it  changes  speed  a  lot  but  stays  on 
the  same  heading. 

For  the  case  of  bearing-only  reports,  however,  the  natural  extension  of 
this  procedure  would  not  give  meaningful  estimates  of  the  driving  noise 
intensity  parameters  because  of  the  less  direct  observability  of  the  motion 
state.  It  would  probably  also  perform  badly  for  what  is  almost  the  same 
thing,  position  reports  with  long,  narrow  "error  ellipses."  The  alternate 
procedure  described  below  can  be  used  in  such  cases  to  estimate  the  driving 
noise  intensity  as  if  the  in-track  and  cross-track  components  were  equal  as 
well  as  independent. 


In  the  notation  of  Refs.  1  and  2,  the  driving  noise  intensity  matrix 
components  and  q  are  now  presumed  equal,  and  are  denoted  simply  by  q 


here . 


The  q  component  is  now  zero, 
xy 


Hence,  only  the  single  scalar  driving 


noise  parameter  q  is  estimated  in  this  alternate  method.  If  the  components  of 

•  R  and  the  orientation  of  the  corresponding 


the  covariance  matrix 


M 


sigma-ellipse  are  as  shown  below. 
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2 

where  M  is  defined  as  in  the  preceding  section,  then  the  quantities  a  ,  b 
sin  g  and  cos  g  can  be  determined  as  follows: 


set 

2 

a  = 

gxx  +  gyy 

set 

ii 

C\J 

X) 

-  e,y)2  *  “h,2 

set 

a2  = 

(a2  +  b2) 

3et 

cos  g 

a2  +  g 
°xx 

,  2  .2 
(a  “  ^  +  8xy 


set  sin  3 


where  "="  is  used  in  the  sense  of  FORTRAN.  Now  let 


which  has  zero  mean  and  covariance  matrix 


M 1  R  +  Qt  . 


The  last  term  in  this  covariance  matrix  is  due  to  an  independent  driving  noise 
increment  with  covariance  matrix 


Qx  !  0  1  „„  „  .  P*  !  0 

.0  T'qx]  •  so  Q  =  [oTq,  ‘ 


Now  let 


=  n 


x  ...  „  r  cos  s  !  Sin  si 

t  L-sin  3  i  cos  6J 


This  rotated  innovation  vector  has  zero  mean  and  covariance  matrix 


a(M1  +  R)nT  +  nQnT  = 


r 

o 

cvj 

CO 

1 - 

r< 

i 

IO 

Lo  !  b2. 

i 

O  1 
l 

< 

T  . 


since  the  contribution  of  the  driving  noise  increment  is  statistically 
independent.  Hence,  if 

2  1  /  2  2, 

s  =  -  (e  -  a  ) 

T  * 

and 

2  1  ,  2  h2. 

n  =  -  (e,  -  b  ) , 

T  ’ll 

these  two  statistics  are  independent  of  each  other  and  their  counterparts  at 
other  report  times,  and 

E(s2)  =  E(n2)  =  q. 


Hence,  the  alternate  procedure  advocated  here  is  to  use  h  =  max  (o,q)  as 
a  running  estimate  of  q  -  together  with  the  other  parts  of  the  tracking 

2  2 

algorithm  as  before  -  where  §  is  the  average  of  the  s  and  n  statistics  for 

2 

all  the  currently  available  reports  (except  that  only  n  is  used  for  the 
bearing-only  reports,  and  the  "e^*"  limiting  procedure  of  Ref.  1  is  used). 
This  average  §  is  computed  recursively  as  before.  Now  there  are  two  recursion 

2  2 

updates  of  §  at  each  instance  of  a  position  report  -  one  for  s^  and  one  for  1^ 


-  but  only  one  for  each  bearing  report.  Of  course,  some  modifications  in  the 

2 

above  procedure  for  computing  p  would  be  needed  for  the  case  of  a  bearing- 
only  report.  One  could  either  replace  it  at  the  outset  with  an  equivalent 
position  report  with  a  very  long  variance  component  along  the  line-of-bearing 


10 


1 


2 

(and  then  not  use  the  computed  s  statistic,  of  course),  or  one  could  compute 

e  as  the  cross-bearing  component  of  the  observed  minus  predicted  position  (as 

^  2 
evaluated  from  the  observed  bearing  at  the  predicted  ship  range)  and  b  as  the 

2  2  2 
cross-bearing  component  of  M  plus  r  o  .  where  r  is  the  predicted  range  and  cr 

is  the  bearing  variance.  With  either  procedure,  a  reasonable  value  of  the 
2 

statistic  n  could  then  be  computed. 
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